Fractal tracer distributions in turbulent field 

theories 



Jonas Lundbek Hansen and Tomas Bohr 

i Center for Chaos and Turbulence Studies 

0^ 



- 

X 



The Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen 0, Denmark 



(D 

Q> ' Abstract 

^ We study the motion of passive tracers in a two-dimensional turbulent velocity field 

■ generated by the Kuramoto-Sivashinsky equation. By varying the direction of the 
. velocity- vector with respect to the field-gradient we can continuously vary the two 
\ Lyapunov exponents for the particle motion and thereby find a regime in which the 

■ particle distribution is a strange attractor. We compare the Lyapunov dimension 
. to the information dimension of actual particle distributions and show that there is 

good agreement with the Kaplan- Yorke conjecture. Similar phenomena have been 

^ ' observed experimentally. 

o 

■ 1 Introduction 



H ' The velocity field of particles confined to the surface of a fluid does not have 

to be incompressible even though the fluid itself is incompressible. It has cor- 
respondingly been observed experimentally [1,2], that tracer particles moving 
only on the two-dimensional surface of a three-dimensional fluid can lie on a 
strange attractor as shown in FIG. 1. In these experiments some 50 million 
floating particles were advected on the surface of the fluid, which was heavily 
stirred from time to time, with the stirrings done so seldom, that the fluid had 
time to come to rest. Pictures of the particle distributions were then exam- 
ined and were found to have a well-defined fractal dimension between 1 and 
2 and hence to be a strange attractor, which changed shape for each stirring. 
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This has been explained theoretically in terms of "random maps" [3,4] , which 
seems appropriate for this pulsed flow. On the other hand, this phenomenon 
is rather general and should also occur in systems with no separation be- 
tween a stirring and relaxing phase. To investigate this possibility we have 
looked at the advection of passive scalars in a simple turbulent field theory, 
the Kuramoto-Sivashinsky (KS) equation, which has been derived e.g. in the 
context of chemical turbulence [5] and the flow of a falling fluid fllm [6] . 

2 Advection by the ID Kuramoto-Sivashinsky Equation 

The Kuramoto-Sivashinsky equation in 1-1-1 dimensions is given as 

hf — hxxxx hxx ~l~ {hx^ (1) 

where /i is a scalar fleld and /j. = The boundary conditions are periodic 
on X e [0, L] . The equation for the derivative u — hx yields 

~ '^xxxx '^xx ~l~ '^UUx- (2) 

This equation has the same nonlinear term as the Navier-Stokes equations. 
It is thus natural to look at the field u = as a. velocity field and study 
the motion of a passive particle with the equations of motion ^ = ahx{x,t), 
where a is a real constant that gives the ratio of the time scales of the particle 
motion and the h-Held. Note that the h-Held itself is not a good choice for 
a velocity field, since the equation (1) is invariant under the transformation 
h ^ h + const. 

Such a study was carried out by Bohr and Pikovsky [7]. They examined the 
case a — 1 and found the particles to diffuse anomalously, i.e., the average 
displacement in time goes approximately as < x{t)'^ >^(x t^/^, which is faster 
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than usual diffusion. They also found, for systems seeded with several parti- 
cles, that particles coalesce in time. Therefore the Lyapunov exponent for the 
particles is negative even though they move in a chaotic velocity field. Thus, 
the system exhibits Eulerian but not Lagrangian chaos. 



3 Advection by the 2D Kuramoto-Sivashinsky Equation 

It is trivial to generalize the Kuramoto-Sivashinsky equation (1) to two di- 
mensions. The equation is given as 

^ = - V^/i - V'/i + ( V/i) • ( V/i) . (3) 

The /i- field is now a function h — h{x, y, t), and V = -|- j^. For the higher 
order derivatives, we use the identities = V ■ V and = V ■ V(V ■ V). 
The reason for the simplicity in the generalization is that all terms are scalar 
products and hence scalars. Again we assume periodic boundary conditions 
on e [0,L] x [0,L]. 

With random initial conditions and periodic boundary conditions, the solu- 
tions look like mountain ranges: mountains (maxima) are separated from each 
other by valleys (minima). The valleys form a web where the mountains are 
isolated holes in the web. In the course of the dynamics mountains are created 
spontaneously or by splitting large mountains into (2,3 or 4) smaller ones - or 
they can disappear by being "squeezed" by neighbouring mountains. FIG. 2 
shows solutions to the 2D KS-equation with L = 50 at four consecutive time 
steps. 
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3. 1 Choice of velocity field 



The most obvious choice of velocity field is, as in the ID field propor- 

tional to the gradients of the /i-field, i.e., 

Vg = oV/i (4) 

where the real constant a, as in the ID-case, is the ratio of the time scales of 
the particle motion and the /i-field. The gradient always points in the direction 
of the local maxima. This causes particles with different initial positions to 
coalesce at these maxima, similar to what is seen in the ID-case [7]. Thus, 
particles moving in the gradient velocity field have both Lyapunov exponents 
Ai, A2 less than zero. The original field h can be regarded as a velocity potential. 
This motion of particles is area contracting and hence dissipative. In terms of 
fluid dynamics, the velocity fleld V/i is compressible. 

We can generalize the equations of motion (4) by a linear function of the 
gradient vectors. Since we want to preserve the isotropy of the KS-equation, 
the only possibility for the generalization is scaling and rotation. The general- 
ized velocity field is thus constructed as a multiplication of a constant a, the 
rotation matrix and the gradient vectors: 



V = aB^tV/i where 'Qq — 



^ cos 9 sin 9 ^ 



sin 9 cos 9 



(5) 



Of special interest is the vector rotated 90° from the gradient: 

Wh = aBgooV/i =1 I . (6) 



This vector field viewed as a velocity field is incompressible, i.e. V • Vh — 0, 
and the scalar field h can be regarded as a stream function. Equivalently, 
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particle motion in this velocity field can be said to be a non-autonomous (time- 
dependent) Hamiltonian system with one degree of freedom, where x and y are 
conjugate variables and h is the Hamiltonian. A Hamiltonian system is area 
preserving and hence the two Lyapunov exponents have the same numerical 
value with opposite signs, (Ai = — A2). Thus, if a volume element is expanding 
in one direction it should contract in the other to preserve the volume [8]. 
Advected in this velocity field, the particles move around the maxima of the 
field at a nearly constant height. The maxima can be regarded as vortices of 
the velocity field. 

Throughout this paper, we restrict ourselves to examining the behavior of 
particles with the most natural ratio of the time scales, thus a — 1. We then 
describe the vector field v by just one parameter, the angle 9. 

3.2 Particle Trajectories 

In FIG. 3 wc show the trajectories of single particles for different values of 9 
advected in the same /i-field. 

It is seen how the trajectory corresponding to pure Hamiltonian motion {9 ~ 
90°) is made up of segments of smoothly rounded curves joined in sharp cor- 
ners. This is because in the changing field, the particle for some time spirals 
around one maximum. When this maximum disappears, the particle changes 
direction drastically and spirals around another maximum. 

The particle moving in a pure a gradient velocity field {9 = 0°) moves at a 
(local) maximum and thus follows the cell's motion. When the cell is destroyed, 
it jumps to a new cell. 

It is seen how the particles travel significantly longer in the Hamiltonian {9 — 
90°) case than in the gradient {9 = 0°) case. 
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3.3 Strange attractor 



For certain values of 9 the particles actually move on a strange attractor. The 
range of the values of 9 that gives a strange attractor can be estimated by 
calculating the Lyapunov exponents for the particles. Given the two Lyapunov 
exponents of a dynamical system, of which at least one should be positive, the 
Lyapunov dimension Dl is given by [8] 

= l + ^ (7) 

The Lyapunov exponents are calculated in the following way [9] : We linearize 
the equations of motion of a small disturbance 5x(t) around the trajectory 
x(t). The equations of motion for the system are given as: 

Xi^Vi{x,t) (8) 



We expand the equations of motion of a small disturbance 5x(t) in a series 
around x(t) to first order and get 

6^n{t) = J(x) • 5x„(t) (9) 

where J(x) = ^g^^t)^ is the Jacobian of the equations of motion in the point 
x(t). We use equations for two disturbances, 5xi(t) and 5x2 (t), since we want 
to find two Lyapunov exponents. The initial conditions for the first distur- 
bance vector are chosen in a random direction with length |5x„(i)| = 1. The 
second disturbance vector is chosen orthonormal to the first. The equations of 
motion of the disturbances are integrated along with the equations of motion 
of the trajectory. If the system is chaotic, at least one of the vectors 5x„(^) 
grow exponentially and it is necessary to re-orthonormalize the vectors from 
time to time to avoid numerical overflow. The re-orthonormalization is done 



6 



in a way such that the first vector is simply renormalized while the other 
vector is turned and renormalized. Let r be the time between consecutive re- 

orthonormalizations. At each re-orthonormalization instant jr, the length aj^' 

(2) 

of the first vector and the area aj spanned by the two vectors are stored. 
The Lyapunov exponents are then given by 



The Kaplan- Yorke conjecture states that Dl actually gives the information 
dimension Di of the attractor of the system. 

We see that for the chaotic Hamiltonian case which has one positive and one 
negative Lyapunov exponent with the same numerical value, — 2, and thus 
the attractor of the particles is a two dimensional subset of the whole plane, 
as expected for a Hamiltonian system. It is thus necessary to introduce some 
dissipation (i.e. 9 < 90°) to find a strange attractor. 

The information dimension is defined as 



where Pi is the natural measure of particles in box i. The natural measure in 
a box is the fraction of the total number of particles, or normalized density, in 
each box. The sum is taken over all boxes i, which has nonzero natural mea- 
sure. In FIG. 4 numerically calculated values of the two Lyapunov exponents 
are shown as a function of the angle 9 (the two lower curves) . It is seen that 
9 has to be quite large, 9 ~ 72°, to obtain one positive Lyapunov exponent. 
Also shown is the Lyapunov dimension of the corresponding attractor. 




(10) 



Di^- lim 




(11) 
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It should be noted, that the Lyapunov exponents and their signs are also 
strongly dependent on the ratio of the time scales, i.e., the values of the con- 
stant a. The Lyapunov exponents decrease with increasing o and the 9 at 
which the transition Ai > occur increases. In the limit of o ^ oo, the transi- 
tion newer occurs although both Lyapunov exponents in the pure Hamiltonian 
system (9 — 90°) are zero corresponding to a h-Held that is constant in time 
i.e. h = h{x, y). 

3.4 Results for the KS Equation 

In FIG. 1 we show pictures of the strange attractor on the surface of a real 
fluid, taken from [1]. It is seen that there are not that many large scale struc- 
tures, the large bends of the strange attractor are big compared to the image 
size. We thus expect the system size of the Kuramoto-Sivashinsky equation 
needed to produce a strange attractor to be only a few cells big. It of course 
has to big enough to be chaotic. A system of size 32 x 32 spatial units were 
chosen. This corresponds to approximately 4 cells in each direction. In the ex- 
periments ~ 50 million particles were used [1] . We thus have to advect many 
particles to see the strange attractor. As a compromise between computing 
time and the large number needed, 190000 particles were used. The parti- 
cles were originally placed in a square grid covering the whole system. The 
2D-Kuramoto-Sivashinsky equation was integrated using an explicit finite dif- 
ference scheme with grid spacings of half a spatial unit, thus 64 x 64 mesh 
point were used. The gradient field between the mesh points was found using 
a bicubic spline interpolation [10]. The computing time for the integration of 
the system with particles was one hour per time unit on the used computers 
(Pentium 166 MHz.) 

After a transient time the particles produce a pattern which resembles a 
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strange attractor which changes in time. Patterns for different values of 9 
but the same /i-field at the same time are shown in FIG. 6-7. 

We have estimated the information dimensions of the patterns numerically, 
using the information dimension Di (11), where the 32 x 32 field is divided 
into N square boxes with side-length e. The dimensions are found as the slope 
of a straight line fit to a plot of Y^f pilogpi versus loge. 

The measured information dimensions with error bars from the fit, the Lya- 
punov exponents and the Lyapunov dimensions calculated from (7) are shown 
in FIG. 4 as functions of 9. 

The Lyapunov exponents change continuously with 9 and with that the infor- 
mation dimension and we can thus choose the dimension of the attractor by 
varying 9. 

In FIG. 5 we show how the measured value of the information dimension 
changes in time. The horizontal lines in the plot are the corresponding cal- 
culated Lyapunov dimensions D^. For the three highest used values of 9, the 
information dimensions are seen to converge nicely to the calculated Lyapunov 
dimensions, while for the lowest used value of 9, the information dimension is 
decreasing in time. This is probably a numerical problem connected with the 
existence of local regions where the velocity field is contracting in all direc- 
tions for a long time, which causes particles to coalesce and thus lowers the 
information dimension. 

4 Conclusion 

We have found fractal tracer distributions of particles advected in the velocity 
field of the 2D-Kuramoto-Sivashinsky equation. The measured dimensions of 
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attractors are in good agreement with calculated Lyapunov dimensions. It 
could be interesting to see if the fractal tracer distributions could be obtained 
for particles advcctcd in other equations exhibiting spatio-temporal chaos, e.g. 
the Complex Ginzburg-Landau equation, where some work on the motion of 
advected particles has been done [11]. 
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Fig. 1. . 

Pictures of tracer distributions from experiments. Left: laminar case. Right: 
fractal case. From Sommerer [2]. 
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Fig. 2. Solutions to the 2D Kuramoto-Sivashinsky equation at four consecutive 
times. Time going with At = 1 from upper left to lower right. The system size is 
50x50. 
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Fig. 3. The trajectory of an advectcd particle for different values of 6 in the same 
/i-field. Prom upper right to lower right = {0° , 18°, 72°, 90°}. The particle moved 
for 90 time units. Prom the numerical solution of the 2D Kuramoto-Sivashinsky 
equation, the field h is known only at discrete points in space. To find the values of 
the velocity field at all points in space, making the /i-field a continuous function we 
used bicubic spline interpolation [10] . This is a standard method for interpolating a 
function which is only known at discrete points. The interpolated function returned 
has continuous derivatives. 
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Fig. 4. The two Lyapunov exponents of the particle motion for different values of 9 
are shown as the two lower curves. The upper curve show the Lyapunov dimension 
as calculated from the Kaplan- Yorke conjecture when one Lyapunov exponent is 
positive. The data points are the measured dimensions of the attractors from the 
simulations 
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Fig. 5. The measured information dimension at different times for each of the four 
used values of 6. The horizontal lines are the corresponding calculated Lyapunov 
dimensions Dl. It is seen how the values of the measured dimensions apparently 
converge to the calculated Lyapunov dimensions for the three highest values of 6. 
For the lowest value of 9, near the threshold of chaos, the dimension of the particle 
distribution apparently does not converge in time. The fits could scldomly be done 
over more than one decade, due to cross-over effects at each end of the interval. 
The limitations are mainly due to computional problems of handling more than the 
190.000 particles used. 
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Fig. 6. A realization of the particle patterns at time t = 60. Values of 9 from upper 
left to lower right: 76.5°, 81.0°, 85.5°, 87.75°. The measured dimensions are from 
upper left to lower right 1.309 ± 0.005, 1.601 ± 0.008, 1.873 ± 0.006, 1.957 ± 0.005. 
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Fig. 7. A realization of the particle patterns at time t = 180. Values of from 
upper left to lower right: 76.5°, 81.0°, 85.5°, 87.75°. The measured dimensions are 
from upper left to lower right 0.979±0.007, 1.49±0.008, 1.821±0.004, 1.942±0.006. 
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